CREATED USING THE RSC LaTeX PCCP ARTICLE TEMPLATE - SEE www.rsc.org/electronicfiles FOR DETAILS 

ARTICLE TYPE www.rsc.org/xxxxxx | XXXXXXXX 

Computational Methodologies and Physical Insights into Electronic 
Energy Transfer in Photosynthetic Light-Harvesting Complexes 



Leonardo A. Pachon"" and Paul Brumer 



*b 



^ I Received Xth March 2012, Accepted Xth XXXXXXXXX 2012 
O First published on the web Xth XXXXXXXXXX 2012 
^ DOI: 10.1039^000000x 
>^. 

^( We examine computational techniques and methodologies currently in use to explore electronic excitation energy transfer in the 
^^ context of light-harvesting complexes in photosynthetic antenna systems, and comment on some new insights into the underlying 
1-H physics. Advantages and pitfalls of these methodologies are discussed, as are some physical insights into the photosynthetic 
CN dynamics. By combining results from molecular modelling of the complexes (structural description) with an effective non- 
^_^ equilibrium statistical description (time evolution), we identify some general features, regardless of the particular distribution in 
,^ |the protein scaffold, that are central to light-harvesting dynamics and, that could ultimately be related to the high efficiency of 
Qn-the overall process. Based on these general common features, some possible new directions in the field are discussed. 
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1 Introduction 

The transfer of electronic excitation energy is ubiquitous in 
nature and its dynamics and manipulation is of special interest 
in diverse fields of physics, chemistry, biology and engineer- 
ing. Traditionally, electronic excitation energy transfer (EET) 
has been studied and described in, e.g., molecular crystals ^ 
rare gases in the solid phase ^'^ and dye aggregates^. More re- 
cently focus has turned to EET in biological pigment-protein 
complexes present in light-harvesting complex of natural pho- 
tosynthetic antenna systems ^ . The dynamics of the EET pro- 
cess in these systems is often thought of in terms of donors 
and acceptors of excitation, and it is characterized, mainly, by 
two time scales: a time scale Tda at which the transfer between 
donors and acceptors takes place (intermolecular transitions) 
and a second time scale Tyr at which the intramolecular (vibra- 
tional) relaxation occurs^. 

Based on the Tda and Tyr time scales, the field distinguishes 
between two limiting cases. If the intermolecular transitions 
'are faster than the intramolecular relaxation, Tyr/Tda < 1 , then 
the exciton is delocalized among the donor- acceptor pairs 
and it travels across the complex coherently. If, by contrast, 
Tvr/Tda > 1, then the exciton is localized and it hops from 
donors to acceptors, and propagates incoherently across the 
complex. Under in-vivo conditions, such as those in photo- 
synthetic light-harvesting complexes^, EET has been largely 
considered as incoherent, and its high efficiency (< 100%) has 
been discussed in terms of the ultrafast, tens of picoseconds, 
character of the overall process ^'^ . [An important exception 
to the incoherent-dynamics assumption is found in, e.g., the 
Light-Harvesting Complex II (LHCII is present in most of the 
higher plants) or in the Chlorosomes (found in, e.g., green sul- 



fur bacteria) where, due to the strong coupling between chro- 
mophores, the excitations stay delocalized over a certain num- 
ber of pigments, even at room temperature ^'^.] 

Recent experimental results ^^~^^ have suggested that dur- 
ing the first steps of EET in some light-harvesting complexes, 
quantum coherences between different monomers of the pho- 
tosynthetic chromophores may be assisting in the high ef- 
ficiency of the electronic excitation energy transfer. This 
provocative suggestion resulted from the observation that 
these coherences were, apparently, living longer (~ 800 fs at 
77 K and ^ 400 fs at 277 K) than could have been predicted 
(cf. Refs. 15,16 or Sec. 5.2 below). These two issues, and 
their implications, rapidly became the center of considerable 
study from diverse fields of chemistry and physics (cf. the 
review articles in Refs. 17-22). 

The initial attempts to explain the dynamics observed 
in these experiments ^^^^ were based on Forster's the- 
Qj,y 23,24 (^stj^ong exciton- vibrational coupling compared with 
the dipole-dipole interaction between the donor and accep- 
tor, fast intramolecular relaxation, incoherent transfer) or on 
the Redfield formalism ^^ (weak exciton- vibrational coupling, 
slow relaxation). These studies revealed that these two ap- 
proaches were inadequate because, in the dynamics of pho- 
tosynthetic complexes, there is neither a dominant time scale 
nor a leading coupling strength ^^~^^. This fact motivated the 
use of approaches and techniques developed in other fields 
such as quantum optics ^^~^^, semiclassical molecular dynam- 
ics^^~^^, solid state physics ^^'^^'^^, etc., and encouraged the 
development of a variety of new techniques ^^'^^ . In addition, 
it motivated the improvement of molecular modelling tech- 
niques ^^'^^'^^, in particular in aspects related to the influence 
of the protein environment, and to the possibility of all-atom 
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simulations of the system time evolution ^^. 

Our goal in this Perspective is to summarize the state-of- 
the-art in the computational description of EET in photosyn- 
thetic complexes as well as to suggest some new possible di- 
rections in the field. In doing so, we collect common features 
that have emerged from the different approaches. We also at- 
tempt to collect the relevant physical information that could 
be used as input for molecular modelling strategies in the po- 
tential design of artificial EET systems. 

2 Photosynthetic Complexes 

Structures of natural light-harvesting antenna complexes are 
extremely diverse. One can find anything from strcutures 
that are highly symmetric structures to those where the chro- 
mophores seem to be placed randomly in the protein scaf- 
^Q^^ 8,3 1,36 jjjjg diversity can be evidenced in, e.g., the in- 
terchromophore separation, which is, e.g., ~ 1 nm for organ- 
isms with chlorophyll-based networks or more than ~ 2 nm 
in organisms with bilin-based networks ^'^^'^^. Remarkably, in 
all cases the light harvesting efficiencies in these networks are 
reported to be close to 100% 831,36 

Despite this variability in structure, the chromophore exci- 
tation energy/electronic coupling landscape share some sim- 
ilarities that are worth analyzing in order to identify mecha- 
nisms that underlie the high efficiency of the photosynthetic 
transport. We postpone this description until Sec. 5 and, in 
the following, describe some of the structural features of two 
of the better studied complexes, the Fenna-Matthews-Olson 
(FMO) present in green sulfur bacteria and the Phycocyanin- 
645 (PC645) complex from Chroomonas sp., in order to pro- 
vide a physical picture of these systems. Subsequently, we de- 
scribe the standard theoretical framework for studying EET ^'^. 
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2.1 Examples: FMO and PC645 

The FMO pigment-protein complex from Chlorobium 
tepidum^^'^^ is a trimer consisting of identical, weakly inter- 
acting monomers ^^. Each weakly interacting FMO monomer 
contains seven coupled bacteriochlorophyll-a (BChl-a) chro- 
mophores arranged asymmetrically, yielding eight nondegen- 
erate, delocalized molecular excited states (excitons)^^'^^. 
The coupling values of the seven chromophores, organized as 
{BChk 1, BChk 2,..., BChk 7}, can be summarized by the 
following Hamiltonian 
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where the energy of the third chromophore has be taken as the 
reference zero. Values are given in units of cm~^ in the site 
representation^^. 

Since the electronic coupling between the chromophores 
BChk 1 and BChla 2 is relatively strong compared to the other 
coupling strengths ^^, see //fmo. it is common to focus, for 
purposes of model studies of EET, on the dynamics of the ex- 
citation in a dimer composed by BChla 1 and BChla 2. (This 
often adopted model is discussed further in Sec. 5). 

Although the seven site model has been the standard, an 
eight bacteriochlorophyll-a chromophore was recently discov- 
gj,^^ 38,39 jj^ what follows we restrict our discussion to the 
case of seven chromophores and refer the interested reader to 
Ref. 33 for a detailed analysis of the influence of the eighth 
chromophore on the EET dynamics. 

The PC645 pigment-protein complex has been studied both 

experimentally ^^ and numerically in great detail ^^ It contains 
eight bilin molecules covalently bound to the protein scaffold. 
A dihydrobiliverdin (DBV) dimer is located at the center of 
the complex and two mesobiliverdin (MBV) molecules lo- 
cated near the protein periphery give rise to the upper half 
of the complex absorption spectrum. Excitation of this dimer 

may initiate the light harvesting process ^^. The coupling val- 
ues of the eight bilin molecules, arranged as {DBVc, DBVd, 
MBVa, MBVb, PCBcl58, PCBdl58, PCBc82, PCBs82}, can 
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be summarized by the following Hamiltonian ^^'^^ 



^PC645 = 
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in units of cm~^ The electronic coupling between the closely 
spaced DBVc and DBVd molecules is ~320 cm~^ This rel- 
atively strong coupling results in the delocalization of the ex- 
citation, yielding the dimer electronic excited states labelled 
DBV+ and DBV— . Excitation energy absorbed by the dimer 
flows to the MBV molecules which are each 23 A from the 
closest DBV, and ultimately to four phycocyanobilins (PCB) 
that absorb in the lower-energy half of the absorption spec- 



trum 
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Note that FMO and PC645 play a distinctly different role 
physically. The former serves to transmit energy from the an- 
tenna to the reaction center. That is, under natural conditions 
light absorption is the function of the outer antenna structure 
and the probability of photon absorption in the FMO is very 
low. By contrast, PC645 does, naturally, play an active role in 
absorbing incident radiation. 

2.2 Hamiltonian 

Theoretical modelling of the exciton energy transfer pro- 
cess comprises intramolecular and intermolecular contribu- 
tions (cf. Chap 9 in Ref. 6, or Ref. 5) and involves a number 
of assumptions, noted below, regarding the presumed physics. 
For the particular case of a pigment-protein complex com- 
posed ofN pigments, the corresponding Hamiltonian, /fppc, 

N 
m=\ m,n>m 

where Hm describes intra-pigment contributions and Vmn con- 
tains all the inter-pigment Coulomb interactions. The Hamil- 
tonian of the m-th pigment is characterized by the electronic 
coordinates r^ and the nuclear coordinates R^ and can be de- 
composed as Hm = Trn + ^m , whcrc Tm denotes the 
nuclear kinetic energy contribution. Since we are not con- 
sidering electron transfer among different pigments, we can 
expand the pigment-protein complex Hamiltonian in terms 
of the adiabatic electronic states {(pmaijm\^m)} defined by 

^m^^(Rm)<Pmfl(rm;Rm) = ^ma9mfl(rm;Rm), whcrC a rcfcrS tO 

the a-\h state of pigment m and Umai^m) denotes the corre- 
sponding single-molecule potential energy surface (PEF). 

The inter-pigment interactions in Eq. (1) depend on the dis- 
tance Yjnn between the electrons in the pigment m and those in 



the pigment n, as well as on distance R^„ between the electron 
in pigment m and the nuclei of pigment n. These interactions 
comprise the electron-electron interaction V^'^^ and electron- 
nuclei interactions V^^„~^"^. As discussed in Sec. 2.2.2.1, the 
coupling to the nuclear degrees of freedom can be identified 
as the first source of electronic decoherence. 

The transfer of electronic excitation energy is often accom- 
panied by the transfer of electrons ^'^^. However, it is gen- 
erally assumed that the pigments are sufficiently far apart 
so that the intermolecular wave functions do not overlap 
{{^maWnb) = ^ma.nb) • Under thesc circumstanccs, the inter- 
molecular interchange can be neglected. As such, we are ne- 
glecting the bridge-mediated electron transfer (cf. Ref. 41 
or Chap. 7 in Ref. 6) that could induce long-range super ex- 
change of electrons ^^ . Nonetheless, one has to be aware that 
strong coupling to the solvent degrees of freedom can induce 
overlapping of the intermolecular wave functions. 

We shall also restrict ourselves to the situation when the 
electronic ground state Sq = g and the first excited singlet 
state 5"! = ^ of the different pigments suffice for describing the 
dynamics of EFT. This is the well known two-level model ^. 
Additionally, any non-adiabatic transitions between 5*0 and S\ 
during FFT are neglected. 

After projecting the Hamiltonian in Fq. (1) onto the 
set of states defined by the Hartree product ansatz 

{nLl 9ma^(rm;Rm)}, WC gCt^ 



^ppc — 2^ L^ ^ma I ^ma ) \ Vma \ 
m a=g,e 

\ 2^ Jmn I ^me ^ng ) \ ^ne ^mg \ i 



(2) 



where the matrix elements H^a = ( (pma |^m| (pma ) and Jmn = 

\ Twe H^ng I ^mn \ Yne Yfng ) • 

From the pigment-protein complex Hamiltonian in Fq. (2), 
it is clear that the electronic energies Hma as well as the exci- 
tonic couplings Jmn depend upon the nuclear coordinates. As 
discussed below, they also depend upon the external environ- 
ment (solvent), so that, in general, their values are constrained 
by the underlying interactions. Methodologies and strategies 
for obtaining these parameters are discussed in Sec. 4. 



2.2.1 The site representation — The transfer of electronic 
excitation energy can be described in terms of the num- 
ber of excitations present in the complex. If no excitations 
are present, then the complex is in its ground state, |0) = 
Hm I Vmg ) • In the case when a single excitation resides at pig- 
ment m, it is denoted \m) = \ cpme ) Hriy^m I %g ) • This ordering 
of states suggests the following decomposition of the Hamil- 
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tonian in Eq. (2), 






= £/f™g|0>(0| 

m 

Y.\^^e^lL^ng\ |m)(m| + ^^„|m)(^| (3) 



where the second Hne in Eq. (3) denotes the Hamiltonian of 
the zero-exciton manifold and the third Hne denotes the single- 
exciton manifold Hamiltonian. In Eq. (3), the intermolecular 
electrostatic coupling has been neglected^. In general, the de- 
scription of the dynamics in terms of the /t-exciton manifolds 
is known as the site representation. 

2.2.2 Exciton- vibrational interactions — As already dis- 
cussed, in addition to the solvent, the vibrational degrees of 
freedom comprise contributions from the single-pigment vi- 
brations and from inter-pigment vibrations ^. The role of these 
vibrational degrees of freedom (DOE) is twofold: apart from 
modulating the electronic energies and the excitonic couplings 
in Eq. (2) or in Eq. (3), they are a source of electronic deco- 
herence in cases where the vibrations are not measured. As 
explained below, in such a case the electronic degrees of free- 
dom constitute the system, whereas the vibrations serve as the 
environment. Note that Sec. 2.2.2.1 contains a fulldiscussion 
of how the nuclear degrees of freedom can induce electronic 
decoherence. 

Note that, following Refs. 5 and 6, we also effectively in- 
clude the effects of the vibrations as well as the influence of 
the solvent (environment or reservoir) around the pigments. 
In Sec. 2.4, we discuss the characterization of this effective 
model of the vibrations and the environment from a statistical 
point of view. 

2.2.2.1 Electronic decoherence by vibrational DOF. As 
an example, consider the excitation from an initial prod- 
uct state comprising the ground electronic state | g ) and the 
ground vibrational state | Vg ) to an excited electronic state | e ) 
in a single pigment. Eor a sufficiently fast laser pulse, the state 
of the pigment will be described by the density matrix 

p,o,=A^(\g)\y,){y,\{g\ + \c\^\e)\v,{t)){v,{t)\{e\ 



A5th 



+^*|g>|v,)(v,(0|(e|e"''« + c|.)|v,(0)(v,Kg|e-" 



\5th 



(4) 



where A is a normalization constant, c is proportional to the 
dipolar element de^g and 5 is the difference in energy between 
the minima in the nuclear potential on the ground and excited 
electronic states. 



Since our interest is in the electronic dynamics then, as- 
suming the Condon approximation is valid, we trace over the 
vibrational states and get the density matrix for the electronic 
degrees of freedom: 



p=A^\g){g\ + \c\^\e){e\+c*\g){v,{t)\v,){e\ 



A8th 



^c e 



){8\{yg\vg{t))c-^'^'). 



(5) 



The electronic coherence, manifested in the off-diagonal ma- 
trix elements \g){e\ and \e){g\, is modulated by the over- 
lap (Vg(^)|Vg). As the vibrational wave packet |Vg(^)) 
evolves on the upper electronic surface, particularly when the 
wavepacket comprises many states, the overlap (Vg(^)|Vg) 
can decrease and so will the off-diagonal terms; the result is 
electronic decoherence (cf. Refs. 42 and 43 or Chap. 6 in Ref. 
44). 

This process, electronic decoherence by the vibrational 
states, has been widely studied using quantum/classical 
methodologies ^^'^^'^^~^^ and it is expected to occur on ul- 
trashort time scales. Eor example, results on betaine dye 
molecules ^^ and on femtosecond dynamics and laser control 
of charge transport in trans-poly acety lene ^^'^^ suggest that 
these time scales are very short, ~2.5 fs and ~3.7 fs, respec- 
tively. We return to this subject in Sec. 5, because it is the 
comparison of these time scales to those observed experimen- 
tally, which are far longer, that has motivated a great deal of 
research. 

In our description we have assumed that the each pigment 
is affected only by its local vibrations, which are assumed to 
be uncorrelated. However, which degrees of freedom are to 
be regarded as the bath and what is the nature of the vibra- 
tions is crucial to the decoherence issue ^^. Eor example, it has 
been shown that if the correlations between local vibrations 
are taken into account, they could suppress the cross-over to 
incoherent dynamics at high temperatures ^^ or affect the trap- 
ping time in rings of chromophores ^^. This subject is widely 
discussed in Ref. 51 and more recently in Ref. 52, references 
to which the reader is referred. 

2.2.2.2 Inclusion of the intra/inter vibrational DOFs and 
the solvent. As a first approximation, one could assume that 
the nuclear motion as well as the low-frequency solvent co- 
ordinates can be treated in the harmonic approximation and 
attempt to simulate its effects via a set of dimensionless vibra- 
tional normal mode coordinates {qa} of frequency (Oa- In this 
case, the Hamiltonian for the ground state contribution of the 
pigment-protien-complex in Eq. (3) now reads 



H, 



^p^ = ^o|0)(0| 



^Q=EQ^Y.^mg^Y.^(Oabiba 



(6) 

(7) 
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where £"0 denotes the ground state energy, including vibra- 
tional zero-point energies, of the normal-mode vibrations^. 
ba and ba are the annihilation and creation operators, re- 
spectively, expressed in terms of the normal mode coordinates 
Qa = {bi^ba) with qa = ^JhJlc^aQa' 

Assuming that the ground and the singly excited electronic 
states can be described by the same normal coordinates, the 
single excitonic Hamiltonian reads ^ 



^w^ = Y^^mn\m){n\ 



m,n 



'^mn — Ofnntifne ~r (1 <^mwj 



^L +X^^a|m^(Of)2o 



L ^kg ^Y.^(Oa[ blba + - + 8m{cc)Qc 



(8) 



(9) 



k^m 



The exciton- vibrational coupling matrix gmn{cc) = 
Sm,ngm{o^) + (1 " Sjn^n)gmn{oc) in Eq. (9) is given by the 
dimensionless coupling constants gm{oc), determined from 
the linearization of the potential energy surfaces along the 
normal modes, and by the coupling constants gmn{cc), which 
account for the influence of the vibrational motions in the 
electronic couplings /^^. In the Frank-Condon approxima- 
tion, vibrational motion does not induce electronic transitions; 

that is, gmn{CC) = Sm,ngm{CC). 



2.3 Summary of the assumptions and approximations 

It is useful to summarize the approximation made in deriv- 
ing Eqs. (7) and (9), which is essentially the Frenkel-exciton 
model. Assumption one is the absence of inter-pigment elec- 
tron transfer, which is guaranteed if there is no intermolecular 
wave function overlap, e.g., if the pigments are sufficiently far 
apart. Second, the EET dynamics is assumed to be fully char- 
acterized by the electronic ground state Sq = g and the first 
excited singlet state 5"! = ^ of each pigment, yielding " the 
two-level model". 

In the next step one couples the electronic manifolds to 
the vibrational nuclear modes and to the solvent coordinates, 
which induce decoherence and dissipation. For simplicity, 
we have assumed no inter-pigment electronic-nuclear interac- 
tions, so that nuclear vibrations pertain to each pigment. Ad- 
ditionally, we assume that the collective nuclear- solvent vibra- 
tions can be described by a global manifold in its harmonic ap- 
proximation. Since the dipole-interaction matrix elements will 
depend, in general, upon this collective coordinate, we further 
adopt the Frank-Condon approximation and neglect any pos- 
sible non-adiabatic transitions induced by the solvent. This set 
of approximation results in Eqs. (7) and (9). 



2.4 Dissipation and Decoherence 

Since the number of nuclear and solvent degrees of freedom 
(DOE) in Eq. (8) is extremely large, one could invoke a sta- 
tistical treatment of those DOE. From a statistical viewpoint 
the combined effects of the vibrational states and the sol- 
vent in each pigment in Eq. (8) can be described in terms 
of a local fluctuating random force Cm(0 describing the fluc- 
tuations of the position-normal-mode coordinate. Each of 
these fluctuating forces can be characterized in terms of its 
two-time correlation function {l^m{t)^m{^))mg, where Qmg de- 
notes statistical average over the equilibrium density operator, 
pmg = exp(— //;^gj3)/tr[exp(— /f^gj3)], of the normal modes in 
Eq. (9). For the particular case described in Eq. (9), the corre- 
lation function is be defined by ^^~^^ 



(Cm(OCm(0))mg = ^ / M^) 

Jo 71 



coth 



dco 

71 



cos(a)^) —i^m{(Ot) 



(10) 



where j3 = \/k^T, T is the temperature of the combined 
environment and Jm{co) is termed the spectral density (see 
Sec. 2.4.1 below). The fact that (Cm(OCm(0))mg is complex 
is a consequence of the anticommutativity of the normal mode 
coordinates pa and qa^^- In the classical limit, H ^ 0, the 
correlation function is real and becomes 



iUtKmm 



mg 



— / —/^ 0) cos 0)0 = -^r^ 0, 

no Jo CO p 



(11) 
where Tmit) is the relaxation function. Tmit) is the directly 
observable quantity in time-dependent fluorescence Stokes 
shift experiments ^^'^^ and is defined in terms of the reorga- 
nization energy A^ as r^(0) = 2hXm^^ . From Tm{t), the 
relaxation time T^ of the environment can be determined as 



r,„(o) 



S^dtT^{t). 



2.4.1 The spectral density. — The spectral density con- 
tains the relevant information about the nature of the envi- 
ronment and is usually chosen to fit some measurable quan- 
tity, e.g., the optical spectrum, of the system. As discussed 
in Sec. 4, spectral densities may also be obtained from com- 
putational methods (cf. Refs. 37,58). Here again the idea is 
that the dynamical features, usually classically described, be 
in agreement with the optical measurable quantities. The com- 
plexity of obtaining reliable spectral densities is indicated by 
the studies in Ref. 58, where a host of highly structured spec- 
tral densities are shown for EMO. In the following, we restrict 
attention to the spectral densities already used for studying the 
dynamics of EET. 

For the particular case of the EMO complex, e.g., Adolphs 
and Renger^^ derived a spectral density which allows calcu- 
lations that show agreement between the experimental optical 
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spectrum and calculations based on pigment transition ener- 
gies. This spectral density is given by 

J¥Mo{co) = co^Sogo{co) + co^Su5{co - coh), (12) 

with ^0 = 0.5, ^H = 0.22, cor = 180 cm-^ The spectral 
density Jfmo{co) comprises two parts, one related to pro- 
tein vibrations and characterized by go{co) = 6.105 x 10~^ x 



(coVcoi4)e-v'^^ + 3.8156 x 10"^ x {coycoz^)^-^^^, 
C0\ = 0.575 cm~^and CO2 = 2 cm~^ and the other associated 
with a vibrational mode characterized by the Dirac delta func- 
tion in Eq. (12). This spectral function was recently used in 
Ref. 34 to study the dynamics of the FMO complex and the 
results are discussed in Sec. 3.2.2. A comparison of this spec- 
tral density function and that computed recently ^^ is provided 
in Sec. 4 below. 

Another type of spectral density widely used in the context 
of EET, as well as in non-linear spectroscopy, is the Ohmic 
model with Lorentz-Drude regularization 15,18,21,56,57,59 



Jmico) = 2hXmTm(0/ (l ■ 



-co^^l] 



(13) 



For this case, the relaxation function reads r^(^) = 
2/zA^exp(— ^/t^). a variant to the Lorentz-Drude regular- 
ization model is the Ohmic spectral density with exponential 

decay ^^'^^'^^ 



Jyyi[(o) = 2/ra)exp(-a)/a)c), 



(14) 



where the dimensionless parameter K describes the damping 
strength and (0^ is the cut-off frequency. For this variant, the 
reorganization energy is given by Xm = IKmhcOc^m and the re- 
laxation time by T^ = ;r/(2ft)c,m). The Ohmic spectral density 
with exponential decay is a member of a more general family 
of spectral densities parametrized by s as J{co) ~ co^e"^/^'^. 
Here, 5- = 1 is the Ohmic spectral density while 5- < 1 (5- > 1) is 
termed the sub-Ohmic (super Ohmic) spectral densities ^^. Al- 
though the Ohmic spectral density is a useful choice for, e.g., 
electron transfer dynamics or biomolecular complexes ^'^^ , it 
has the drawback that it does not contain any high-frequency 
modes of the environment. ^4,58, 62 

Another relevant spectral density is the one characteriz- 
ing the effect of blackbody radiation (e.g., natural sunlight or 
moonlight) ^^~^^ . From an open-quantum-system perspective, 
the influence of the blackbody radiation is condensed in the 
spectral density ^^~^^ 



Jbb{co) = Mtbb co^ ^Ib/ (^bb + ^ 



(15) 



where Me = m + MeTBB^BB is the renormalized mass of 
the electron (whose bare mass is m), Tbb = 2e^/3MQC^ ~ 
6.24 X 10~^^s and ^bb is a frequency cutoff ^^~^^. For this 
spectral density, the decay function in Eq. (11) is FbbI'^') = 
Tbb^bb P^('^) ~ ^BBexp(— ^bbI'5'I)] . In the limit ^bb -^ °°, 



one gets the surprising result Fbb (s) = 0. This corresponds 
to taking the point-electron limit. An interesting feature of 
Eq. (15) is that it is capable of introducing fluctuations with- 
out dissipation, with no violation of the fluctuation-dissipation 
theorem ^^. As such, it provides a sound statistical mechani- 
cal formulation for the pure-dephasing approaches (cf. Refs. 
27,28). We have explored the excitation of open quantum sys- 
tems by blackbody radiation in Ref. 68 and some associated 
comments are given in Sec. 5. 

3 Spectral-Density Based Approaches 

The fact that in photo synthetic light-harvesting complexes 
there is neither a dominant time scale nor a leading cou- 
pling strength ^^~^^ means that modelling the excitation energy 
transfer demands the use of robust techniques beyond con- 
ventional theoretical approaches such as Forster's theory ^^'^^, 
or Redfield^^ or Lindblad equations ^^'^^ . Due to the vari- 
ous approximations employed in these methodologies, such 
as the use of perturbation theory in the system-environment 
coupling strength ^^'^^'^^, the incoherent transfer approxima- 
^JQj^ 18,77,78^ the assumption og Markovian dynamics ^^'^^, or 
the secular approximation ^^'^^, these approaches can provide 
inaccurate descriptions of light-harvesting dynamics. They are 
therefore of limited applicability when treating realistically 
parametrized models of light harvesting complexes ^^. 

Before discussing recent approaches, it is worth discussing 
the most basic theory for EET: the Forster theory ^^'^^. As al- 
ready noted above, the theory developed by Forster applies 
to the case of incoherent EET, i.e., to the case when the inter- 
molecular transfer times are slower than the intramolecular re- 
laxation. This approach is based on the Fermi golden rule and 
on a second order approximation to the excitonic coupling be- 
tween pigments. This implies that Forster's theory describes 
EET for weak exciton coupling. Forster's theory characterizes 
EET in terms of excitation rates km^n from pigment m to n. 
These rates are given in term of the spectral overlap between 
the donor emission and acceptor absorption spectra and deter- 
mined by ^'^^'^^ 
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where Fm{co) = J^^d^e^^^e-^^^--^^-^-^-^') denotes the 
emission spectrum of the pigment m and An{co) 
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-gn(t) denotes the absorption spectrum of the 



pigment n. Here, gm{t) is the line-broadening function 18'57,79 
Within the second-order cumulant expansikon (not required in 
the original Forster theory), gm{t) is given by: 
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This expression clearly shows that the dynamics, as well as 
the transport properties, depend on the statistics of the bath as 
well as on the non-local (in time) character of the correlation 
function. 

We note that at present, in the related open-quantum-system 
community, there is a great interest in quantifying the non- 
Markovian (non-local in time) character of dynamics of given 
systems. In doing so, some measures for the degree of non- 
Markovian behavior have been given in Refs. 80 and 81. An 
interesting connection betweemn our discussion above and 
that subject can be made as follows. The intermediate integral 
over ^2 in Eq. (17) accounts for the non-Markovian character 
of the dynamics. This characteristic feature could be useful in 
order to experimentally quantify the non-Markovian character 
of a particular system based on its line-broadening function. 
Work along this direction is currently in progress ^^. 

Below, we divide the discussion into two main parts. In 
Sec. 3.1, we discuss some of the main methodologies based 
on master equation descriptions while in Sec. 3.2 we discussed 
those based on the propagating function. 

3. 1 Master-equations 

Master equation approaches provide equations for the prop- 
agation of the density matrix for the system dynamics (here 
the electrons), obtained by tracing over the environmental 
degrees of freedom. Popular master-equations based ap- 
proaches for studying EET are the second-order perturbative 
time-convolution (TC2) and the second-order perturbative - 
time convolutionless (TCL2) quantum master equations (see 
Chap. 9 in Ref. 83 for a formal basis of these techniques) 
and the hierarchical second-order cumulant expansion ap- 
proach ^^'^^'^^. These approaches are based on the projection 
operator technique ^^-^^ , The interested reader may find useful 
the discussions in Refs. 83,85-88. 

A comprehensive and detailed discussion regarding the ap- 
plicability of these three methodologies using typical val- 
ues present in photo synthetic complexes, Jmn ~ 100 cm~\ 
Xmn "^ 100 cm~\ Ejn ~ 100 cm~^ and T^ ~ 50 cm~^ at 77 and 
300 K for the dimer composed of \(pie)\ 92g ) and | (pig )\(p2e)^ 
was presented in Ref. 18. For completeness, we briefly dis- 
cuss only the main features of these approaches. 

- TC2 is calculated in second-order perturbative approxi- 
mation with respect to the coupling of the system to the envi- 
ronment in the interaction picture ^^'^^'^^'^^. Neglecting initial 
correlation between the pigments and the solvent/vibrational 
DOFs, in this approach, the time evolution of the system den- 
sity matrix is calculated as (cf. Chap. 9 in Ref. 83) 

where J^mit^s) is the convolution or memory kernel ^^. The 



equation (18) is a simplified version of the Nakajima-Zwanzig 



equation 



85,86 



Although, the non-local character of this 



integro-differential equation accounts for the non-Markovian 
character of the dynamics, it also makes Eq. (18) cumbersome 
to implement. Thus, further approximations/assumptions are 
needed. 

Using a classical version of the fluctuation-dissipation the- 
orem [equivalent to taking the high temperature limit for the 
real part of the correlation function of the noise in Eq. (10)], it 
was shown in Ref. 18 that the transition rates predicted by the 
TC2 equation deviate strongly from that given by the Marko- 
vian Redfield equation as well as from those predicted Forster 
theory ^^. The authors in Ref. 18 conclude that the TC2 equa- 
tion is applicable only for the nearly Markovian regime, in 
spite of its non-Markovian nature. 

In Ref. 89 it was shown that the Markovian approximation 
to the TC2 gives reliable results only for small reorganization 
energies (A^ < 10 cm~^ for typical values of EET in photo- 
synthetic complexes '-«^ 100 cm~^). This is far from the domain 
relevant to light harvesting. 

- TCL2 is based on the TC2 plus the application of the 
time-convolutionless projector operator technique of Shibata 
et al ^^, which allows for a time local description of the mas- 
ter equation in Eq. (18). TCL2 is accurate for describing 
coherence between electronic ground and excited states in a 
monomer, regardless of the magnitude of the electron-phono n 
coupling ^^. However, as pointed out in Ref. 18, TCL2 fails 
to describe the transfer rate in a region of large reorganization 
energy, A^ > 20 cm~^ For small reorganization energies, the 
rate predicted by the TCL2 equation is virtually the same as 
that for the Markovian Redfield equation ^^. 

- Second Order Cumulant and Hierarchy Expansion. This 
strategy allows for the inclusion of site-dependent reorgani- 
zation energies as well as an appropriate description of the 
Gaussian fluctuations of the bath^^'^^'^^. On the basis of the 
high temperature approximation, this approach has been used 
as a benchmark and as a reference in the field because it is 
able to extrapolate between the Redfield and Forster theories. 
However, it is restricted to bilinear system-bath coupling, and 
places very large demands on memory ^^. Further, the numer- 
ical effort grows rapidly with increasing system size, and with 
spectral densities that deviate from the Lorentz-Drude form^^ . 
Notwithstanding, this approach has the advantage that can be 
implemented to be numerically exact and is currently used by 
many research groups (cf. Refs. 18,92). 

Note that 77 K corresponds to ~ 50 cm~^ and 300 K to 
~ 200 cm~^ These are the typical temperatures at which 
the 2DPE experiments have been done. One can immediately 
see that the high temperature condition fiEm/k^T <C 1 adopted 
above in some approaches above is barely fulfilled for light- 
harvesting photosynthetic complexes {fiEm/k^T ~ 2 at 77 K 
and fiEm/k^T ~ 0.5 at 300 K ). In order to be applicable to 
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other physical scenarios, where the deviations from the high 
temperature are larger, it would be of interest to efficiently ex- 
tend these methodologies to the low temperature regime. 

3.2 Propagating-function methods 

The second main category of approaches relies directly on 
the propagating function of the system and bath density 
matrix. We discuss below three of the most commonly 
used: linearized-density-matrix-dynamics approaches, the 
quasi-adiabatic propagator path integral (QUAPI) and the non- 
interacting-blip approximation (NIB A). 

3.2.1 Linearized-density-matrix-dy namics — Lin- 
earized approaches ^^~^^'^^~^^ for quantum evolution were 
developed in order to deal with non-adiabatic chemical 
processes such as proton and electron transfer in solution, 
excited- state molecular fragmentation or molecular relaxation 
after electronic photo excitation ^^'^^ . In general, these 
approaches assume a quantum system represented by a 
number of discrete basis states | n ) coupled to an, in principle 
highly dimensional, environment described by continuous 
coordinates (Q^P) (cf. Ref. 29). The Hamiltonian for such a 
system, in the diabatic representation, is given by 



The density matrix is evolved according to 



Pit) 
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+ i:/iAv(e)(|A)(A'| + |A')(A|). 

X<X' 



(19) 



The first term represents the nuclear kinetic term while the 
remaining terms comprise the electronic Hamiltonian. Note 
that the Hamiltonian in Eq. (19) is analogous to that in Eq. (8). 
An important ingredient in these approaches 29-32,93-98 j^ ^^ 
map the discrete quantum states onto continuous coordinates. 
The strategy is based on the mapping formalism ^^~^^. The 
idea is to replace the evolution of the electronic subsystem 
with the evolution of a system of fictitious harmonic oscilla- 
tors^^. This is achieved in two steps: (i) Map the electronic 
Hilbert space spanned by the n diabatic states into one coin- 
ciding with a subspace of n harmonic oscillators of unit mass 
and at most one quantum of excitation in a single specific os- 
cillator^^, i.e., I a) -^ \ma) = |0i,... la...,0^), and (ii) re- 
place the projection operator by harmonic oscillator creation 
and annihilation harmonic operators, | A )( A' | -^ a^^d^f, and 
express these in terms of their oscillator coordinates and mo- 
menta {qa^Pa), <^A' = (^A' +i^AO/V^- Once the mapping is 
performed, the original Hamiltonian in Eq. (19) becomes 



^n. = Py2M^Whu{Q) [ql ^Pl-h) 

^ X 



(20) 



Thus, the resulting Hamiltonian is expressed completely in 
terms of operators with continuous spectra. Despite this fact, 
it is able to mimic the effects of transition among discrete elec- 
tronic states ^^'^^. 

The second main ingredient in the different linearized ap- 
proaches is the linearization itself. To do this, the unitary time- 
evolution operator q-'^^^I^ is represented in terms of discrete 
phase space path integrals in the environmental variables and 
double sums over the quantum states ^^'^^. From the semiclas- 
sical perspective, the evolution of the density matrix is gov- 
erned by forward and backward trajectories associated with 
the unitary time-evolution operator and its adjoint operator, 
respectively. The way in which the action is linearized in the 
different approaches is what distinguishes between them; the 
particularities of each are discussed below. 

In the linearized approximation to the semiclassical initial 
value representation (LSC-IVR) of the unitary time-evolution 
operator (cf. Ref. 30), the action in the path integral expres- 
sion is expanded to linear order in the distance between the 
forward and the backward trajectory of the environmental and 
electronic degrees of freedom. Thus, the the effect of the bath 
(nuclear degrees of freedom) on the dynamics of the system 
and the dynamics of the systems itself are approximated. This 
approach was used in the past^^~^^ and more recently used for 
the particular case of the FMO complex by Tao and Miller ^^. 
Despite the level of approximation, their results for the dimer 
case are in agreement with those of Ishizaki and Fleming ^^'^^ . 
However, when the full complex is considered, the results at 
this level of approximation are not reliable ^^'^^. 

In the linearized approach to nonadiabatic dynamics using 
the mapping formalism (LANDmap)^^'^^'^^, the linearization 
is performed only in the environmental degrees of freedom. 
LANDmap has been used for studying the dynamics of the 
FMO complex in Ref. 29 as well as the dynamics of the PC645 
complex in Ref. 3 1 . In particular, results from Ref. 29 show 
that when the bath relaxation is slow (e.g., Tc ^ 500 fs) re- 
sults obtained with LANDmap are reliable for times around 
1 ps at high temperature (~ 300 K) and for a wide range of 
reorganization energies [A ~ — 5A, being A the strength of 
the electric coupling (A = 100 cm~^ in Ref. 29)]. For fast 
relaxation (e.g. Tc ~ 100 fs) and low temperature (~ 77 K) 
LANDmap approach provides a reasonable representation of 
the population oscillations only for short times {r^ 100 fs) and 
for small reorganization energies (A < A/5); for this set of 
parameters, LANDmap is expected to fail to reproduce ther- 
malization^^'^^. 

LANDmap can be used to generate a short time propagator, 
which can be iterated to generate more reliable results in the 
long time regime and valid for a wider spectrum of parame- 
ters; this scheme is known as the iterative linearized density 
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matrix (ILDM) approach ^^'^^. Although ILDM is more accu- 
rate than LANDmap, its convergence after many iterations can 
be problematic ^^ and it is more demanding in computational 
terms. 

The use of the mapping formalism brings its own draw- 
backs, in particular when implemented and accompanied by 
semiclassical approximations (cf. Ref. 32 for details). To cir- 
cumvent these problems, Huo and Coker have developed an 
improved version of ILDM: the partial linearized density ma- 
trix (PLDM)^^. It is based on a coherent state representation 
of the electronic part and a linearization of the nuclear (envi- 
ronmental) degrees of freedom. This approach was tested for 
nonadiabatic multi states scattering problems and it compared 
very well with standard results ^^ . An important feature of this 
new approach is the good performance at long times and the 
correct equilibration with the thermal bath. However, being 
based on a coherent state representation, PLDM suffers from 
an excess of free parameters ^^~^^^ However, the additional 
parameters can be used to improve the convergence of the av- 
erage results. 

The advantage of the linearized approaches is that they can 
be applied, in principle, to any molecular model with arbitrary 
system bath interactions. In particular, PLDM offers a favor- 
able balance between accuracy and computation efficiency ^^. 

3.2.2 QUAPI-based methods — The quasi-adiabatic 
propagator path integral (QUAPI) ^^^~^^^ has been exten- 
sively used by Nalbach and Thorwart in the context of 
gg J 34,105,106 ^g |j^ previous approaches, QUAPI assumes 
that the density matrix of the donor- acceptor system plus 
the environment can be characterized by a product state at 
t = 0, i.e. p(0) = Pda ^Penv The time evolution of Pda(0 i^ 
obtained after tracing out the environmental (or bath) degrees 
of freedom, i.e., Pda(0 = Tr{^(^,0)p(0)t/-H^O)}g , 
U{t,0) being the propagator of the full system plus bath. 
The bath is modelled by a collection of harmonic modes 
^B = 5 L^ {pj + ^i^j) i^ thermal equilibrium at temperature 
T. 

QUAPI is based on a symmetric Trotter splitting of the 
short- time propagator J^{tk^i^tk) for the full Hamiltonian 
into a part depending on the system Hamiltonian and a part 
involving the bath and the coupling term^^'^^^'^^^. The split- 
ting of propagator in terms of the short- time propagator is, by 
construction, exact in the limit 5t = tj^^i — t^ ^ 0. However, 
for finite 5^ it introduces a finite Trotter error, which has to be 

eliminated by choosing 5^ small enough until convergence is 

achieved 34^102,103^ 

Another key ingredient in QUAPI' s performance is the 
treatment of non-Markovian time evolution generated from 
the non-local correlations of the bath modes (see Sec. 2.4). For 
any finite temperature, these correlations decay exponentially 
fast at asymptotic times, thereby setting the associated mem- 



ory time scale 34,102,103 j^^^ f^^^ allows for the introduction 
of an effective memory-time window Tmem = ^St, so it is as- 
sumed that the system dynamics have memory over K slices of 
time. Over this time window, QUAPI defines an object called 
the reduced density tensor, which has to be iterated in order 
to propagate the reduced density matrix of the system. Within 
the memory time window, all correlations are included exactly 
over the finite memory time Tmem and can, in principle, be ne- 
glected for times beyond Tmem- Then, the memory parameter 
K has to be increased, until convergence is found. So, we have 
the requirement of decreasing 5^ while increasing K in order to 
get convergent results: the two strategies naturally run counter 
to one another. However, convergent results can be obtained 
over a wide range of parameter regimes ^4, 102,1 03 

The efficiency and implementation of the QUAPI algorithm 
are based on the choice of the parameters k and M (the num- 
ber of basis states), e.g., the reduced density tensor is a com- 
plex array of size M^^^^. This limits its applicability to the 
case of model systems that are not too large. For the case 
of the donor-acceptor dimer system, M = 2, so with the stan- 
dard hardware architectures one can choose k< 12 — 14^ i'^^^. 
However, if one is interested in studying, e.g., the network of 
eight bacteriochlorophyll-(2 chromophores in the FMO com- 
plex within the two level system framework one has M = 16, 
so K* ~ 3. This relatively short memory window could be prob- 
lematic if one is interested in using this approach in the con- 
text of electronic transfer, e.g., in semiconducting carbon nan- 
otubes^o^. 

Although QUAPI is an iterative scheme, the truncation of 
the memory makes the simulations scale linearly with the 
propagation time; however, since at low temperatures the bath 
induced correlations decay only algebraically ^4, 1 08, 1 09^ ^^^^ 

truncation limits the applicability of QUAPI to the finite tem- 
perature scenario ^4, 102, 103 j^^ comparative studies ^4, 1 06^ |^ j^^^ 

been shown that QUAPI is able to generate the correct time 
evolution of photo synthetic complexes in the range of param- 
eters of interest, finite temperature, strong system-bath inter- 
action and strong electronic coupling. 

QUAPI was implemented by Nalbach et al in Ref. 34 
in order to include the spectral density given in Eq. (12) de- 
rived in Ref. 37. The result suggests that the use of this spec- 
tral density provides slightly smaller coherences lifetimes than 
those observed in the experiment ^\ those calculated using the 
Ohmic spectral density in Eq. (13), and those based on a hy- 
brid quantum/classical all atom calculation^^. It would be in- 
teresting to analyze this case using other methodologies in or- 
der to understand the role of the vibrational high frequency 
modes (see Sec. 2.4.1) in the dynamics. 

3.2.3 NIBA-based methods — In a series of papers, Cao 
et a/. ^^'^^^'^^^ have explored the efficiency of the energy 
transfer in light-harvesting systems using approaches such 
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the Haken-Strobl model (pure dephasing)^^^, the general- 
ized Bloch-Redfield equation (second-order cumulant expan- 
sion model ^^^'^^^) and the non-interacting-blip approxima- 
tion (NIBA)^^'^^. The Haken-Strobl model has been exten- 
sively discussed in the literature ^^ and the generalized Bloch- 
Redfield approach shares the spirit of the Nakajima-Zwanzig 
equation discussed in Sec. 3.1, so we discuss here the NIBA 
approximation ^^'^^ only. 

To understand the term non-interacting blip, we can appeal 
to the path integral formulation of the systems with a finite 
Hilbert space ^^'^^. For simplicity we restrict the discussion 
to the case of a system with two levels. Time evolution of 
double-sided objects, like the density matrix, is performed 
by a pair of trajectories ^^'^^^ In the case of the two level 
systems, these trajectories are piecewise constant paths with 
sudden jumps between states, so we have four possible "po- 
sitions", namely both trajectories visit the ground or excited 
state or one of them visits the ground state while the other one 
is the excited state and vice- versa. The two paths can be seen 
as a single path visiting the four configurations. A period of 
the path spent in the diagonal configuration (ground-ground 
or excited-excited) is called a sojourn while a period the path 
spends in an off-diagonal configuration is a blip^^'^^ . The 
NIBA assumption is that the average time the system spends 
in a diagonal configuration is much larger than time it spends 
in an off-diagonal one. This assumption motivates setting (i) 
the sojourn-blip and (ii) the blip-blip correlations to zero; this 
leads to the non-interacting blip approximation. In the strict 
Ohmic case, i.e., in the scaling limit (Oq -^ oo, assumption (i) 
is exact ^^; thus, the approximation enters, in this case, in ne- 
glecting the interblip interactions. 

An alternate, apparently less complicated, approach to de- 
riving the NIBA is given in Ref. 114. 

According to Ref. 56, this approximation can be justified 
for: 

1 . Weak-coupling to the environment and zero bias. 

2. Super-Ohmic spectral densities with ^ > 1 and ^ > 2 at 
zero and finite temperature, respectively. 

3. Sub-Ohmic spectral densities with ^ < 1 and ^ < 2 at zero 
and finite temperature, respectively. Therefore, NIBA is 
justified ^^'^^ (i) at all temperature in the sub-Ohmic case, 
(ii) at high temperature in the Ohmic case, and (iii)at high 
temperatures in the super-Ohmic case with s < 2.^^'^^. 

NIBA can be formulated as a second-order master equation 
in the bath-dressed electronic coupling ^^^, which is more con- 
venient for the extension to multilevel systems ^^ . This strat- 
egy was followed in Ref. 33 to generate a version of NIBA for 
multi-state systems. This extended NIBA was compared to 
the generalized Bloch-Redfield approach by Cao et ah ^^ for 
describing the time evolution of site populations in the FMO 



complex with eight bacteriochlorophy 11-^2 chromophores; the 
results are in excellent qualitative agreement. NIBA's perfor- 
mance was also studied in Ref. 116 in a wide range of pa- 
rameters. In the limits described above, it provides accurate 
results. 



4 Molecular Modelling of Photosynthetic Com- 
plexes 

Molecular modelling of photosynthetic complexes is of pri- 
mary importance insofar as information about electronic cou- 
plings, site energies, spectral densities, linear absorption spec- 
trum, etc. are obtained from such studies ^^'^^'^^'^^^'^^^. Re- 
cently, Mennucci and Curutchet^^ presented a complete per- 
spective on these approaches and Konig and Neugebauer^^ 
reviewed the advantages and pitfalls of most of these theoret- 
ical methods in providing EET coupling constants. Here we 
review some of the most recent progress and physical insights 
from the molecular-modelling studies. 

As described in Ref. 20, among the various proposals for 
molecular modelling one can distinguish two methodologies 
for describing the environment and its effects on screening the 
Coulomb interaction between the donor and acceptor transi- 
tion dipoles. The first method considers the environment as 
a dielectric continuum. Its effects are characterized by means 
of a screening factor \/n^, where n is the solvent refractive 
index (Forster-like theories). The second approach treats the 
environment at the atomic level, either by using molecular me- 
chanics force fields or by incorporating a full quantum me- 
chanical description of the chromophore-environment. 

Mennucci et al have developed a method, reliant upon 
the dielectric-continuum-based approach, using a combina- 
tion of the quantum linear response and a structureless- 
polarizable-continuum-media model of the environment ^^'^ ^^ . 
The method is able to deal with the non-equilibrium response 
of the system and of the environment during fast processes, 
such as those involved in electronic transitions and electronic 
energy transfer ^^. It has been successfully applied to exam- 
ine the screening induced by the environment in the electronic 
couplings for a set of over 100 chromophore pairs, includ- 
ing chlorophylls, bilins and carotenoids, taken from structural 
models of photosynthetic pigmentprotein complexes ^^' ^ ^^ . An 
interesting outcome from this methodology is the fact that 
photosynthetic light-harvesting can be tuned by heterogeneous 
polarizable environments of the proteins ^^^ and that the final 
resonance energy transfer step could occur on a timescale of 
15 ps. According to Ref. 120, such a rapid final energy trans- 
fer step cannot be reproduced by calculations based on the 
spectral density description, which predict the energy transfer 
times to be on the order of 40 ps^^. 

The second method, where one carries out a quantum all- 
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atoms calculation for the structure and dynamics of photo- 
synthetic complexes, is currently beyond computational reach. 
In order to overcome this problem, some methods based on 
quantum/classical ^^'^^^~^^^ and continuum/atomistic descrip- 
tions have been proposed ^^'^^^. For a review of the main diffi- 
culties in the implementation of these approaches, such as the 
inclusion of the polarizability of the solvent and the hetero- 
geneous character of the environment, the interested reader is 
referred to Refs. 20,22. 

Apart from the difficulties in computing the various elec- 
tric couplings using atomistic an description, a major chal- 
lenge is the full-atom description of the dynamics itself. A 
first step towards an all-atom calculation of the dynamics of 
photosynthetic complexes, in particular of the FMO complex, 
was presented by Aspuru-Guzik et al. in Ref. 35. Based on 
the Bom-Oppenheimer approximation and adopting the Con- 
don approximation, they decomposed the total system Hamil- 
tonian into three parts: a system Hamiltonian acting on the ex- 
citon sector (described by a set of two-level systems), a bath 
of vibrational modes, and a coupling term between them; this 
corresponds to the discussion in Sec. 2.2.2.2. The molecu- 
lar energies are computed using time-dependent density func- 
tional theory with the dynamics done by the Wigner method, 
i.e. the exciton dynamics are done quantum mechanically and 
the bath dynamics are done classically. The evolution of the 
excitonic density matrix is obtained from a statistical ensem- 
ble of unitary evolutions obtained by solving a time-dependent 
Schrodinger equation. This approach allows consideration of 
nonidentical fluctuations across all the sites. The population 
dynamics of the chromophore is in accordance with previous 
calculations ^^'^^'^^. Based on this approach, Aspuru-Guzik et 
al. ^^ concluded that the site energy cross-correlation between 
chromophores does not play a significant role in the energy 
transfer dynamics. This is in accord with the results of 01- 
brich in Refs. 123,124 for the FMO complex and with those 
of Huo and Coker^^ for the PC645 complex. 

It is worth mentioning that the statistical-ensemble- 
averaging approach has inspired the application of quantum 
state diffusion approaches such as those discussed in Refs. 
126,127. In particular, it has been used to study the role of 
quantum oscillations as well as the dependence on site ener- 
gies in electronic excitation transfer in the FMO complex ^^^, 
the influence of the vibrational modes ^^^ and the influence of 
noise, disorder, and temperature on localization in excitonic 
systems ^^^. These techniques appear accurate, efficient and 
are valid for a wide spectrum of parameters 1^8-130 

In addition, we note that direct connections are being es- 
tablished between all-atoms calculations and master equation 
approaches that use spectral densities ^^'^^^ . For example, Ref. 
58 obtained spectral densities associated with FMO by such 
a full-atom computation. Results, whose effect on FMO dy- 
namics is yet to be determined, show structured /(ft)) with 



considerably higher values of ft) contributing to the spectral 
density than assumed in the Ohmic or Drude-Lorentz models 
or Adolphs-Renger model ^^ previously utilized. 

Finally, note that recent evidence has been found for corre- 
lated fluctuations of site energies and intersite electronic cou- 
plings as well as electronic-electronic coupling that could be 
more significant than the apparently uncorrected site energy 
fluctuations ^^^'^^^ Hence, an interesting possible direction in 
this regard is to develop technique which provide spectral den- 
sities that include bath-dependent electronic terms ^^^. 

5 Insights into the Observed Coherences 

Below we discuss some of the issues relating to the role of the 
observed EFT dynamics in natural photosynthetic processes, 
and introduce an analytic model that provides physical insight 
into the origins of the observed long-lived coherences. 

5.1 Natural processes 

We remark on two issues in need of considerable clarification 
if we are to understand the role of the observed coherences 
in nature. First, note that under natural conditions, photosyn- 
thetic complexes are excited by sunlight, an incoherent source 
of light ^^~^^ that is incident on the system for times that are 
huge compared to the time scale of the molecular dynamics. 
As a consequence, excitation of this kind cannot generate co- 
herent dynamics among the pigments, neither in the case of the 
unitary dynamics associated with isolated molecules, 63,64,66,67 
nor in the case of open systems, where the molecule is in con- 
tact with an environment ^^'^^. Based on this argument, it is 
not yet clear whether coherences induced by coherent fem- 
tosecond laser pulses can play a role under in vivo conditions. 
Second, there remains a concern that the observed coher- 
ences may not necessarily be solely electronic in nature, i.e., 
suggestions have been made that vibrational coherences are 
being observed as well (cf. Ref. 133). Only very recently 
has an experimental protocol been applied that, based on the 
nature of the observed 2D Photon Echo experiment, can dis- 
tinguish between electronic and vibrational coherences. A re- 
cent application of this protocol ^^^ for the case of PC645 has 
identified one of the observed long-lived coherences as elec- 
tronic, and the others as likely vibrational. In this regard we 
note that vibrational coherences would be expected to deco- 
here slowly ^^^ so that understanding such a coherence feature 
presents no significant qualitative challenge. 

5.2 The origin of long-lived coherences 

The major part of this paper has provided a review of a num- 
ber of techniques and approaches for the study of the EFT 
dynamics of photosynthetic complexes. Examination of these 
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detailed computations shows, however, that they do not pro- 
vide direct insight into the basic physical origin of the ob- 
served longevity of the coherences. Below we address the 
question as to the physical origin of the long-lived observed 
coherences, by introducing an analytically soluble model that 
provides correct results for dimers that characterize the FMO 
and PC645 systems. 

At high temperatures and for weak coupling to the environ- 
ment, a classical treatment of the thermal fluctuations with an 
Ohmic spectral density with exponential decay (see Sec. 2.4.1) 
predicts that the electronic coherence discussed in Sec. 2.2.2.1 
decays at a Gaussian rate that can be determined as ^^'^^^'^^^ 



Tg = \h'/2Xk^T, 



(21) 



where X is the system reorganization energy. Based on this 
expression, the dephasing time for photo synthetic complexes, 
for a typical value of A = 130 cm~^), can be estimated ^^ to be 
Tg = 45 fs at r = 77 K and Tg = 23 fs and T = 294 K. 

By contrast, the experiments in photo synthetic com- 
plexes such as the FMO complex ^^'^"^ and the PC645 com- 
plex ^^, discussed above, have found that electronic coher- 
ences among different chromophores survive up to 800 fs 
at 77 K^^ and up 400 fs at room temperature^^. Gener- 
ally longer coherence times have been observed for FMO in 
Ref. 14. This surprising observation and its possible con- 
sequences for biological processes have been discussed ex- 
tensivelyii'i^-i^'2^'2^'31'^4,35,84,110,138 ^^^ ^iw^^ provided mo- 
tivation for the development of the methodologies discussed 
above. Interestingly, despite the diversity of approaches and 
techniques, most^^'^^'^^'^^'^^'^^'^^'^^'^^^ predict long-lived co- 
herences on the same time scales as those found experimen- 
tally ^^'^^'^^. This suggests that the underlying physical fea- 
tures are correctly contained in these approaches. However, 
the sheer complexity of these computations has limited one 
from identifying these essential physical features ^^. 

In order to explore the physical features responsible for 
the survival of these coherences, we (in Ref. 16) discuss the 
case of the relatively strongly coupled dimer composed of the 
BChlla 1 and BChlla 2 in the FMO complex (see Sec. 2.1), 
and the dimer formed of chromophores DVBc and DVBd in 
the PC645 complex (see also Sec. 2.1). The generic Hamil- 
tonian for the two two-level system is described in terms of 
Pauli spin matrices by^^'^^ 



h h 

+ -5^110-^,1^1 + -5a120,,2^2 +^1 +^2, 



(22) 



where Ri = Y,a ^oc,i ( <^a,/ + ^a n ^^ ^^^ reaction field operator 
for molecule /, Bi = Y,a ^^a,/^a i^oc.i is the energy stored in 




Fig. 1 Left hand side: The pair of interacting chromophores. Right 
hand side: the effective light harvesting two-level system formed 
from the pair of interacting chromophores 61416 



the solvent cage of molecule / and 5jij is the difference be- 
tween the dipole moment of the chromophore j in the ground 
and excited states ^^'^^ . The first two terms in Eq. (22) are the 
contributions from the individual sites and the third term is 
the A coupling between them. The subsequent terms describe 
the system-bath coupling. Following Refs. 61 and 116, the 
Hamiltonian in Eq. (22) can be written with respect to the ba- 
sis {\gi) 1^2), \8\) ^ \ei), \e\) (g) 1^2), \e\) ® Vi)} describing 
the state of the two chromophores, i.e. 

i=l,2 a 

I -(£++y+) \ 

h -(e_ + y_) 2A 

^2 2A e_+y_ 

V e++V+ ) 

where e± = ei ± £2, and y± = ?>\i\R\ ± ?>\i2Ri. 

Since under excitation by weak light only the singly excited 
states need to be taken into account, we can identify ^^'^^^ 
the active environment coupled 2D-subspace as {|^i) 
1^2), 1^1) ® k2)}- In this central subspace of Eq. (23), the 
effective interacting biomolecular two-level system Hamilto- 
nian reads 






h 
2' 



a 
i=\2 



(23) 



where e = e_ and V = y_ . This is schematically illustrated in 
Fig. 1, where A is the associated "tunneling energy", between 
the new basis states \e\) \g2) and \g\) |^2)- 

Following the description given in Sec. 2.2, we as- 
sume that the two baths are uncorrected {R\{t")R2{t')) = 
{R2{t^')Ri {t')) = 0. Hence, Eq, (23) can be written in the stan- 
dard form of the spin-boson model ^^ 



H 



[he 



fi 



M(jJ+-a,^g^(/?^+/?J)+^to^Z?U^ 



(24) 
where the h^ includes harmonic oscillators coupled to both 
chromophores, with couplings g^ . The effects of the baths are 
treated in terms of the Ohmic spectral density with exponential 
decay presented in Sec. 2.4.1. 
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The Hamiltonian in Eq. (24) differs from the one described 
in Sec. 2.2.2.2 by the fact that in Eq. (24), both states \ei)^ 
I g2 ) and \g\) 0\e2) are coupled to the same effective bath 
modes, whereas the treatment presented in Sec. 2.2.2.2 shows 
that for the case of a pair of donor-acceptor pigments, the 
ground and excited states should be coupled to independent 
modes. Although, this fact introduces fluctuations of the refer- 
ence energy of the donor- acceptor system, if one is interested 
only in the dynamics of the Bloch oscillation this is irrelevant 
(cf. Ref. 5 1,59). For strong coupling to the bath or high degree 
of non-Markovian dynamics, one expects the two models, e.g, 
are driven to different thermal states. However, the fact that 
our results, see below, coincide with the overwhelming major- 
ity of the current results in literature, means that the coupling 
strength and the degree of non-Markovian character, for the 
particular cases discussed below, are still in a regime where 
both models coincide. Note also that the previous predictions 
made for the lifetimes (cf. Ref. 15) were based on the spin 
boson-model. Thus, in order to get an idea about the lifetimes 
of the coherences, the model is quite well justified. 

The parameter range within which the dimers described 
above, B Chi la 1 and B Chi la 2 in the FMO complex and the 
chromophores DVBc and DVBd in the PC645 complex, lie 
allows for the use of the non-Markovian NIBA plus first or- 
der corrections in the interblip correlation strength, i.e. an en- 
hanced NIBA approximation. For a description of the range of 
applicability of the bare NIBA cf. Sec. 3.2.3. The enhanced 
approximation is valid for weak system-bath coupling and for 
8/2A < 1, over the whole range of temperatures (see Chap. 21 
in Ref. 56), and provides simple and accurate analytic expres- 
sions for relaxation and decoherence rates. As it was discussed 



K 


£/2A 


kjiT/2hA 


2A/fflc 


hT/hOc 


e/coc 


0.105 


0.428 


0.305 


1.052 


0.321 


0.45 


0.105 


0.428 


1.098 


1.052 


1.154 


0.45 



Table 1 Parameters used for dimer formed of BChla 1 and BChla 2 

at r = 77 K (first row) and T = 277 K (second row). 

in Sec. 3.2.3, in the limit when (O^ is much larger than the other 
frequencies in the system, NIBA is exact. In the present case, 
see Table 1, this condition is not fulfilled in all cases; however, 
we do have here an improved approximation which provides 
lifetime estimates that are in very good agreement with exper- 
imental as well more-refined-theoretical results. 

The high temperature limit in this approach is given by tem- 
peratures well in excess of 7b = ^(A^^ + e^)^/^//:B, where 



where in our case A = 2A. Here the spectral density is ohmic 
[Eq. (14)] where K describes the coupling strength to the bath 



and cOc is the bath cutoff frequency. 

For the set of parameters listed in Table 1 , we find 7b ^ 
288 K. Hence the FMO experiments, at 77 K and 277 K, are 
in the low temperature regime, T < T\y. In this regime, the 
Rabi frequency Q,, the relaxation rate jr and the decoherence 
rate 7 are given by ^^ 



^' 






-2KAlf^ [9ti/A(iMb/27r^B^) - ln(Mb/2;r^B7^)] 
= 7r7rcoth(Mb/2^B7^)Aeff/Ab, 
--r,/2^2nK{£^/Al)kBT/h, 



(25) 

(26) 

(27) 



respectively, where Ab = a/A^^ + e^ and \i/{z) is the digamma 

function. (Note that in Refs. 139 and 56, the term A^^^ is miss- 
ing in the expression for Q..) The expressions in Eqs. (25)- 
(27) were derived for the particular case of non-Markovian 
Ohmic dissipation, however, analogous expressions for arbi- 
trary spectral densities can be found in Refs. 16,56,139. 

Based on this model, and using the parameters in Table 1 
(obtained from experimental results ^^), we find^^ the fol- 
lowing time scales for the FMO complex 27tQr^ = 163 fs, 
y-^ = 90 fs, 7-1 = 153 fs at r = 77 K, while 27:0.'^ = 151 fs, 
Y~^ =45 fs, 7~^ = 69 fs at r = 277 K. For the case the dimer 
formed of chromophores DVBc and DVBd in the PC645 com- 
plex, we predicted 2nQ.~^ = 49 fs, 77^ = 76 fs and 7"^ = 88 fs 
at r = 294 K using the parameters in Refs. 16,29,40. Note 
that Y~^ and 7"^ are related to decay rates and that the time 
period over which coherence is observed is therefore much 
longer ^^. For example, electronic coherences are observed to 
persist in PC645 at room temperature for times up to 400 fs, 
a result predicted by this model as well. Similarly, the time 
scales over which coherences persist for FMO, predicted from 
this model, are in agreement with all other (far more complex) 
computations, and with the earlier experimental results in Ref. 
11. 

Hence, the long time scales seen to emerge naturally here 
from the system parameters. Why then were far shorter de- 
coherence time scales originally expected for these systems? 
To see this, note that in molecular systems, dynamics is often 
studied between different electronic eigenstates of the system, 
separated by greater than ~ 10^ cm~\ with no coupling be- 
tween them. In such cases, the dephasing time from Eqs. (26) 
and Eq. (27) would be extremely short. By contrast, in the 
case of photosynthetic complexes, energy transfer occurs be- 
tween exciton states that are close in energy and, additionally, 
are coupled. This generates a small value for the ratio e/2A 
which in turn is responsible for longer dephasing times (see 
Fig. 2). 

Additionally, expressions such as Tq [Eq. (21)], which are 
often used to estimate rates, are only valid at high temperature, 
k^T ^ hcOc, and at short times, t < co~^. Under conditions 
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Fig. 2 Decoherence time 7~^ (continuous blue line) as functions of 
e/2A, for fixed A = 35 cm"^ and at T = 77 K. The dashed blue line 
denotes the dephasing rate used in Refs. 14,28,137 given by 
70 = 2n(k^T /h)X/h(Oc while the dotdashed blue line denotes the 

decoherence time Tg = Jfi^ /iXk^T used in references 15,136,137. 

The relaxation time y^^ is depicted by the black dotted line. Fixed 
parameters as in Table 1 . The data for this plot was taken from Ref . 
16. Note that color is online, but not in the printed text. 



rally and are not surprisingly long. This physical picture de- 
scribed above, was subsequently verified by Cao et ah ^^ using 
the Hamiltonian described in Sec. 2.2.2.2. A similar physical 
picture for the survival of coherences has emerged indepen- 
dently in other physical systems such as spin charge qubits in 
quantum dots^^^'^^^ or in electronic excitations in semicon- 



ducting carbon nanotubes 
to the findings in Ref. 16. 



107 



This provides some robustness 



5.3 Coherences and energy transfer 

The possibility that the coherences observed in the 2D elec- 
tronic spectroscopy experiments ^^^^ could be assisting the 
high efficiency of the photo synthetic process has been one of 
the more exciting conjectures in the field. However, there is 
now sufficient theoretical evidence from different groups, such 
as Whaley et al i^^, Olbrich et aO^^^^\ Cao et aO^^^^^^^^\ 
Coker et al^^ and Aspuru-Guzik et al^^, suggest that the 
presence of these coherences has only a slight effect on en- 
ergy transfer. 



when the expression for Tq is valid, the bath modes can be 
treated classically ^^^, as in Refs. 45 and 46. When one is not 
in the appropriate regime, the classical evolution of the bath 
underestimates quantum coherence effects ^^^ because at low 
temperatures quantum fluctuations overcome thermal fluctua- 
tions ^^^. Hence, estimates based on Tq in Eq. (21) are unre- 
liable. Moreover, as we already mentioned, Tq characterizes 
a Gaussian decay, while in the experimental results i^'i^'^^ i\^q 
decay is observed to be exponential. 

Similarly, the expression y^ = 27z{k^T /fi)X /fico^, another 
usual expression for estimating the exponential decay at high 
temperature, also provides an inadequate representation of the 
true physics and associated dependence on system and bath 
parameters. This expression is derived in the limit £^0,K^ 
and 0)c -^ ^, so it is expected to result in decoherence times 
that are underestimated and misleading ^^. 

In Fig. 2, we have shown the explicit dependence of the var- 
ious decay rates in Eqs. (25)-(27) for fixed values of the ratio 
e/2A using the FMO parameters at 77 K. We observe that the 
decay rate increases for large values of this ratio. Based on 
this result and our previous discussion, one can identify three 
physical features found to be responsible for long coherence 
lifetimes: (i) the small energy gap between excitonic states 
(ii) the small ratio of the energy gap to the coupling between 
excitonic states, and (iii) the fact that the molecular character- 
istics place the system in an effective low temperature regime, 
even at ambient conditions. 

In this framework, the observed long lifetimes arise natu- 



6 Some Future Directions 

Results discussed above apply to the specific case of EET in 
photosynthetic light harvesting. We remark on a number of 
possible related extensions of interest. 

1. It is worth noting a critical and interesting overview 
on EET modelling in light-harvesting complexes by Huo and 
Coker in Ref. 3 1 . They point out that despite the complexity, 
diversity and variability of the protein-pigment complexes (see 
Sec. 2.1), all of these complexes are highly efficient and the 
Hamiltonians of these systems, parametrized using the molec- 
ular modelling methodologies discussed in Sec. 4, share some 
generic features ^^: (i) Clusters of chromophores with closely 
spaced excitation energies that have appreciable electronic 
couplings between the cluster members, (ii) Chromophore 
states whose excitation energies are isolated, but which ex- 
hibit appreciable electronic coupling to neighboring states. 
(iii) These isolated but coupled states are often arranged en- 
ergetically in cascade or barrier patterns that channel the di- 
rectional flow of energy through these multichromophore net- 
works, and toward reaction centers. 

These features, and those discussed in Sec. 5.2, could be 
used as the target properties for artificially designed light- 
harvesting complexes. These characteristics could be the de- 
sirable outcome of molecular designs based on current or new 
molecular modelling strategies, such as those described in 
Sec. 4. Work in this direction is worth considering, mainly 
because it would relate the studies above to the ultimate goal 
of designing artificial light-harvesting devices inspired by na- 
ture, and currently of great interest (e.g., Ref. 144) . 
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2. The physical insights discussed in Sec. 5.2 describe 
how coherence can be preserved after a single photon exci- 
tation. Recently, motivated by possible natural scenarios such 
as photo-protection^, there is interest in the possibility that the 
observed coherence can be restored continuously in time ^^^ . 
This situation suggests a role for non-equilibrium phenomena 
in EET. 

The importance of non-equilibrium contributions can be ex- 
amined in terms of the results derived, in a different context, 
by Pachon et al in Ref. 146. They have shown that for out- 
of-equilibrium quantum systems, the border between the clas- 
sical and the quantum realms is more intricate than that for 
the equilibrium situation. In particular, they have shown that 
some quantum features, such as entanglement, can survive un- 
der higher temperatures in non-equilibrium cases than in the 
equilibrium case^^^ (see also Ref. 147). However, even if the 
process described in Ref. 146 could protect the coherences, 
it is not clear whether or not they could optimize EET, even 
though preliminary results show that these coherence could be 
spatially directing energy transfer ^^^. 

3. The results described above show that the notion of a 
"high temperature regime" is relative, it depends on the en- 
ergetics of the complex and as well as on the solvent prop- 
erties. That is, room temperature could be low temperature 
for some systems but high temperature for others. The fact 
that EET in photo synthetic complexes relies on the moder- 
ate/low temperature can be seen to be a way of protecting the 
possible coherent dynamics that could take place in the light 
harvesting process. The reason for this is that, at low temper- 
atures and in the long time regime, the decay of the correla- 
tions is expected to be algebraic, 1/^^, rather than exponential, 
exp(— 7^^),^"^'^^^'^^^. Thus, in the low temperature regime the 
coherences are expected to decay more slowly than in the high 
temperature regime. Our discussion in Sec. 5.2 offers the basis 
for exploring the design or the search of natural complexes in 
an effective low temperature, where the electronic coherences 
could live even longer than seen at present. 

4. Although methodologies such as QUAPI^"^'^^^'^^^ or the 
second order cumulant expansion ^^'^^'^^ have provided valu- 
able insights into EET, we call attention to some other promis- 
ing new methodologies based on their reliability over a wide 
spectrum of parameters, their performance, and efficiency. 
They are also versatile, applying to different kinds of coupling 
terms and spectral densities. Specifically, we are referring to 
methodologies such PLDM (cf. Sec. 3.2.1) and those ^^^"^^^ 
based on quantum state diffusion ^^^'^^^ that are worth con- 
sidering beyond the context of EET in photo synthetic com- 
plexes, e.g., in superconducting carbon nanotubes ^^^ or other 
physical solid-state physics systems ^~^. After completing this 
Perspective, we noted the current interest in applying and 
developing efficient methodologies based on density-matrix 
renormalization-group approaches (cf. Ref 148,149), which 



seem to be highly efficient and could be used to explore pos- 
sible mechanisms assisting EET in regions of the parameter 
space beyond the scope of our analysis (cf. Ref. 150) or with 
highly structured spectral densities. 
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